Transcript Characteristics on the Susceptibility Difference of Bovine Respiratory Disease

Bovine respiratory disease (BRD) is one of the major health issues in the cattle industry, resulting in significant financial crises globally. There is currently no good treatment, and cattle are made resistant to pneumonia through disease-resistant breeding. The serial blood samples from six Xinjiang brown (XJB) calves were collected for the RNA sequencing (RNA-seq). The obtained six samples were grouped into two groups, in each group as infected with BRD and healthy calves, respectively. In our study, the differential expression mRNAs were detected by using RNA-seq and constructed a protein–protein interaction (PPI) network related to the immunity in cattle. The key genes were identified by protein interaction network analysis, and the results from RNA-seq were verified using reverse transcription-quantitative polymerase chain reaction (RT-qPCR). A total of 488 differentially expressed (DE) mRNAs were identified. Importantly, the enrichment analysis of these identified DEGs classified them as mainly enriched in the regulation and immune response processes. The 16 hub genes were found to be related to immune pathways categorized by PPIs analysis. Results revealed that many hub genes were related to the immune response to respiratory disease. These results will provide the basis for a better understanding of the molecular mechanism of bovine resistance to BRD.


Introduction
Bovine respiratory disease (BRD) is one of the most serious health issues affecting the cattle sector globally, resulting in significant financial losses [1]. In the United States, the incidence of BRD is 70-80% with a mortality rate of 40-50% [2][3][4][5]. In China, the incidence of BRD is 50-100%, and the mortality rate is 20-50%. Extreme weather conditions, such as cold winters and temperature fluctuations between spring and autumn in Xinjiang China, contribute to increased susceptibility of calves to respiratory diseases including BRD [6]. The respiratory illness affects 40-80% of people, with a 60-90% fatality rate. In XinJiang, BRD is caused predominantly by infection with bacteria [6]. Unfortunately, the incidence of BRD and mortality rate increases annually despite the heavy use of antibiotics in the cattle industry [7,8].
BRD is typically caused by bacteria and viruses. There are currently no effective vaccines, and the use of antibiotic treatment is also not ideal. Antibiotics can help prevent sickness; however, their overuse can result in the resistivity of these broad-spectrum antibiotics in human pathogens such as methicillin-resistant Staphylococcus aureus [9][10][11]. The selection of calves that are resistant to respiratory illness genetically may be an efficient method to minimize disease incidence [12]. Increased resistance may lead to reduction in the use of broad-spectrum antibiotics and the improved sustainability and appeal of resulting food products. Using transcriptome (RNA) sequencing, we sought to select cattle that were not susceptible to respiratory disease to reduce the prevalence and influence of disease in cattle.
RNA sequencing (RNA-seq) is a powerful tool that can be used to identify differentially expressed genes (DEGs) as potentially important indicators in breeding; analyze the category, structure, and expression changes in the transcription factors; and reveal the molecular regulatory mechanism of specific biological processes. Transcriptome sequencing has been widely used in disease research [13,14]. In the present study, we analyzed blood samples obtained from calves with different level of resistance to BRD using RNA-Seq to identify key genes implicated in immune response. Our results provide reasonable foundation for elucidation of molecular mechanism relevant to resistance of Xinjiang brown cattle calves to respiratory disease.

Materials and Methods
2.1. Animal Assembly and Sample Retrieval. All experiments were performed in accordance with the guidelines established by the Animal Care Committee of Xinjiang Agricultural University (No. 2018017). Three diseased Xinjiang brown calves and three healthy Xinjiang brown calves were purchased from Altay, Xinjiang (Xinjiang, China). The calves (3 weeks old) were grown in equivalent environments, with ad libitum food and water, shared-housing arrangements. Blood samples were extracted from the jugular vein, flash-frozen in liquid nitrogen, and kept at −80°C. The susceptible and resistant groups were then housed separately.
In Xinjiang, China, the main factor that causes bovine pneumonia is Pasteurella multocida [6]. Animals were identified as BRD-positive based on clinical examination and serum haptoglobin concentration. Calves with at least one visual BRD sign, a rectal temperature ≥40°C, and abnormal  International Journal of Genomics lung sounds detected at auscultation were defined as BRD cases. Calves with no visual signs of BRD, a rectal temperature <40°C, no abnormal lung sounds detected at auscultation, and a serum haptoglobin concentration <0.25 g/L were classified as non-BRD animals for transcriptome analysis. At the same time, nasal swabs of healthy cattle and sick cattle were collected for isolation and identification of pathogenic bacteria. Through polymerase chain reaction (PCR) technology and sequencing, it was found that the selected positive sample was infected by the P. multocida type-A, and no bacteria were isolated from the healthy cattle. No drugs were used in the area/farm from which the animals were purchased for the present study.

RNA Isolation and Preparation.
Initially, total RNA from whole blood was retrieved using trizol reagent (Invitrogen Co. Ltd, San Diego, USA) according to the reported kit procedure. However, the Agilent Bioanalyzer 2100 (Agilent Technologies, CA, USA) and the NanoPhotometer® spectrophotometer (IMPLEN, Westlake Village, CA, USA) were used to calculate the RNA purity in pricked samples.
Consequently, after passing the quality test, RNase-free DNase set was further used for RNA purification, the RNA was purified using RNA purification kit (Tiangen, China). Additionally, the purified RNA quantification was assessed through Agilent Bioanalyzer Quality testing. For the following process, total RNA fulfilling the integrity criteria was selected, i.e., RINs >8.5 integrity and a 28S/18S ratio of >0.7.

Preparing Libraries, Assembling, and Sequencing
Transcriptomes. According to the reported protocol, six RNA libraries for input data were generated using NEBNext® Ultra™ RNA Library Prep Kit for Illumina® (NEB, Ipswich, MA, USA). Poly-T oligo-attached to magnetic beads were used to extract mRNA from the prepared RNA library. In the NEBNext First Strand Synthesis Reaction Buffer, fragmentation was performed five times with divalent cations at a high temperature. M-MuLV Reverse Transcriptase and random hexamer primers were used to make first-strand cDNA (RNase H-). Following that, DNA polymerase I and RNase H were used to synthesize second-strand cDNA. The library fragments were purified using the AMPure XP technology to select cDNA fragments ranging from 150 to 200 bp (Beckman Coulter, Beverly, USA). Before PCR, 3 μl of USER Enzyme (NEB, Ipswich, MA, USA) was used with sizeselected, adaptor-ligated cDNA at 37°C for 15 minutes and then 5 minutes at 95°C. Phusion High-Fidelity DNA Polymerase, universal PCR primers, and index (X) primer were used in the PCR. According to the manufacturer's instructions, the prepared library quality was assessed using the Agilent Bioanalyzer 2100 system, and PCR products were purified (AMPure XP system). The TruSeq PE Cluster Kit v3-cBot-HS (Illumina) was used to cluster the index-coded samples on a cBot Cluster Generation System.

Data
Analysis. Accordingly, the six cDNA libraries were sequenced through the Illumina NovaSeq 6000 sequencer. FastQC platform was utilized to evaluate the quality of obtained raw data. However, the original sequences filtered through FastQC were washed and screened by applying the selection criteria of low-quality data, N-containing reads, and linker sequences through Seqtk software. The filtered clean reads were compared to the Bos taurus reference genome using the Hisat 2.0 software (ARSUCD1.2). The expression level of genes was normalized to calculate fragments per kilobase of exon per million fragments mapped (FPKM) value using the program edge R. The resulting P values were adjusted using the Benjamini and Hochberg's approach for controlling the false discovery rate. The findings were screened differently and counted using the P -value <0.05, and log fold change |log 2 FC| > 1 was used as thresholds for significant differential expression in this study. The heat map of RNA expression was produced using the correlation of normalized gene-level FPKM values across samples with the heat map.

Functional Enrichment Analysis of DEGs. The Gene
Ontology (GO) terms enrichment analysis for identified DEGs was performed through GOseq R [15]. The GO enrichment study was carried out by collecting all of the  3 International Journal of Genomics GO keywords significantly enriched, and further screening the DEGs based on their biological activities. All DEGs were mapped to GO words in the database (http://www.geneontology.org/), and gene numbers for each term were determined using the hypergeometric test to produce highly enriched GO terms. DEGs were considered significantly enriched by GO keywords with adjusted P-values less than 0.05. To perform pathway enrichment analysis and assess the statistical enrichment of DEGs in Kyoto Encyclopedia of Genes and Genomes (KEGG) (http://www.genome.jp/kegg/), the KO-Based Annotation System (KOBAS) program was used [15][16][17][18]. This approach was used to find genes involved in metabolic or signaling pathways that were highly overrepresented. The transcriptional factors (TFs) in the cattle genome were identified and classified using AnimalTFDB (http://www.bioguo .org/AnimalTFDB/). 2.6. Core Hub Protein-Protein Identification. The String database (http://string-db.org/) is a tool for determining how proteins interact. We used the String database to import the differential genes and evaluate the protein interaction    Table 1 lists the forward (F) and reverse (R) primer pairs used to amplify the genes of interest in RT-qPCR procedures.
The following conditions were used to carry out the amplification reactions: 95°C for 10 minutes, followed by 40 cycles of 95°C for 15 seconds, and 60°C for 1 minute. The following circumstances were used to conduct the melting curve Gene ratio: the ratio of the number of differential genes annotated to the GO number to the total number of differential genes. BP: biological process; MF: molecular function; CC: cellular component. 6 International Journal of Genomics analysis: to guarantee that a single product was amplified in each reaction, 1 minute at 95°C, 2 minutes at 65°C, and a gradual rise from 65°C to 95°C have been used. The relative gene expression level was evaluated by the 2 −△△Ct method.

RNA Sequenced Data Analysis.
We analyzed the blood samples of susceptible (S1, S2, S3) and resistant groups (R1, R2, R3) using RNA-seq. After stringent quality assessment (removal of connector and low-quality sequences) of raw RNA-seq data, we obtained 309,071,922 and 298,963,120 raw and clean reads, respectively, in six samples. The average percentage of clean reads was 95.98%. The contents of A and T were similar in the sequencing data, and the curves were coincident. The results of C and G were similar to those for A and T. Eventually, the obtained data were aligned against the reference genome. The assembly of these genomes showed the sequenced data were evenly distributed on the reference

Evaluation of DEGs.
Significantly, a total annotated genes present in susceptible and resistant groups consisted of 14,247 and 13,956 genes, respectively. The maximum values of FPKM in two groups were 1:42 × 10 4 and 0:93 × 10 4 , and the minimum values were 3:61 × 10 −3 and 3:06 × 10 −3 , respectively. The above results indicate that the resistant group expressed fewer genes than the susceptible group, and the magnitude of gene expression was also lower (average and peak) in the resistant group.
A total of 488 genes were found according to P < 0:05 and |log 2 FoldChange| > 1 when transcriptome data of two groups were compared. Subsequently, it was observed by RNA-seq data analysis that about 488 genes were differentially expressed (DE) in susceptible and resistant groups. From these, 353 (72.33%) DEGs showed upregulation while 135 (27.66%) DEGs were down-regulated, respectively. All obtained DEGs are highlighted in the volcano plot along with their hierarchical clustering (heat map) as shown in Figures 1 and 2, respectively.

Functional and Biological Processes Enrichment Analysis.
Additionally, the enrichment analysis for identified DEGs was studied through GO terms, classifying these DEGs based on their cellular component (CC), biological process (BP), and molecular function (MF) GO terms. By using cut off value as P < 0:05, we recognized 68 genes significantly upregulated in molecular transducer activity, receptor activity, serine-type endopeptidase, and peptidase activities, carboxylic acid transmembrane transport, organic acid transmembrane transport, and mainly in the immune response-regulating signaling pathway, regulation of immune response as shown in Figure 3 and Table 2.
Gene ratio: the ratio of the number of differential genes annotated to the GO number to the total number of differential genes. BP: biological process; MF: molecular function; CC: cellular component.
Furthermore, metabolic pathway analysis was performed for these identified DEGS utilizing the KEGG server. It showed that DEGs were significantly enriched with 52 numerous metabolic and signaling pathways, i.e., fructose and mannose metabolism, cytokine-cytokine receptor interaction pathway, nitrogen metabolism, including the IL-17 signaling pathway, nucleotide-binding oligomerization domain (NOD)-like receptor signaling pathway, influenza A, Tumor Necrosis Factor (TNF) signaling pathway, Nuclear factor-kappa B (NF-κB) signaling pathway, Cyclic adenosine 3 ′ ,5 ′ -monophosphate (cAMP) signaling pathway, and Janus kinase/signal transducers and activators of transcription (JAK/STAT) signaling pathway ( Figure 4 and Table 3).

Gene Interactions Network
Analysis. STRING (http:// string-db.org/) tools were used to predict protein interactions among the 488 DEGs. The circular arrangement in Cytoscape software (http://www.cytoscape.org/; version 3.8.0) was used to distribute all genes ( Figure 5).
Among these DEGs, the top 17 hub genes were selected based on their properties. The network and interaction analysis of these hub proteins were further generated and visualized through CytoHubba plugin of Cytoscape ( Figure 6).

Discussion
The survival ratio of calves has a direct bearing on overall economics of cattle production. Respiratory disease has the highest morbidity and fatality rate in calves [19]. Although there are a few studies [19,20] on the molecular regulatory mechanism of BRD using omics techniques, however, studies on the calf model are still rare. The pathogenesis of BRD has a complex regulatory mechanism. And different pathogens may have various regulatory mechanisms. Our current understanding of its mechanistic basis remains limited. In the present study, we focused on the calf respiratory disease caused by P. multocida.
Pasteurella was previously found to be the main infectious pathogen among calves in the Xinjiang region of China [6]. It is one of the most common pathogenic bacteria associated with BRD potentially causing inflation related to platelet aggregation [21].
We found that most of the DEGs were up-regulated in calves with respiratory disease and were associated with innate immunity. The DEGs may be used to predict respiratory diseases caused by bacterial pathogens and as molecular markers for breeding. The identification of key upstream regulatory factors may provide possible targets for developing new molecular therapies [22]. Sixteen key candidate genes related to immune response and function were found  International Journal of Genomics by RNA-seq analysis of susceptible and resistant groups. Additionally, TLR4, CD59, CD36, NFκBIA, and C5AR1 were revealed as the key regulator genes in BRD through proteinprotein interaction network analysis.
Previous studies have shown the inflammatory response of endotoxin-induced bovine mastitis. Furthermore, it is reported that the NF-κB signaling pathway activated by TFs is mainly involved in microbial infection in mice [23][24][25][26][27][28][29]. These pathways may be activated by the toll-like receptors (TLRs) mainly function in adaptive immunity through T-lymphocytes [30][31][32][33][34][35]. When bacteria invade the lungs of calves, biological signals induce the formation of heterodimers of TLR4 and TLR6, which rapidly trigger an inflammatory response [36]. It is believed that specific immune factors can be produced to protect against pathogens that invade the lungs by inducing the NF-κB signaling pathway [30,32]. In light of the foregoing information, it's worth noting that TLR4 was up-regulated in the current study, implying that TLR4's immunological mechanism was activated in the vulnerable group and may have led to an inflammatory response. In addition, NFκBIA was also identified in this study, which is an important regulatory factor in the NF-κB signaling pathway [31].
SMPDL3B, was one of the 16 candidate genes that were up-regulated in the present study. It indicates that SMPDL3B may have inhibited the function of TLR4 and blocked the immune response in the susceptible group [35]. This may have resulted in the immunosuppression of the susceptible group.
The identified marker genes were observed to be functionally enriched in immune regulation and serine activity classified through GO terms. Inflammation is closely related to immunity, and infection evokes the body's immune response, leading to the stimulation of the host cell immune signaling pathway, promoting the inflammatory response, and increasing the release of IL-6 pro-inflammatory factors [32]. Previous studies have found that after infection with bacterial pneumonia, the concentrations of IL-2 and TNF-α in susceptible types were increased significantly and gradually decreased as inflammation declined [32]. The levels of immune factors may be related to the severity of the disease [26,33]. Increasing evidence has shown activation of immune response in human host by the expression of inflammatory cytokines [34].
KEGG pathway enrichment analysis revealed that the NF-κB signaling pathways were significantly enriched. TLR4 was discovered to be a critical candidate gene that has been reported to activate NF-B through Lipopolysaccharide (LPS) and promote IL-1, IL-6, and IL-8 inflammatory factors overexpression [36]. These factors have a role in innate immune responses and disease resistance, as well as activating innate and adaptive immunological responses and aiding in pathogen protection [26,34].

Conclusion
In the present study, we used RNA-seq to perform a comparative transcriptomic analysis of peripheral blood in Xinjiang Brown Cattle calves susceptible and resistant to BRD. We established the gene expression profile of Xinjiang brown cattle. A total of 488 DEGs were screened out. From these, mRNA expression of 16 key candidate genes related to immune response, i.e., SMPDL3B, IL5RA, NFE2, CD59, CX3CR1, CLEC4E, TLR4, CD36, NFκBIA, C5AR1 GZMB, HGFAC, IL15, TMPRSS2, ALOX15, and CCL24 was validated using RT-qPCR. The RT-qPCR gene expression trend was consistent with the results of RNA-seq. According to core hub, genes were analyzed by degree methods, we infer that TLR4, SMPDL3B, and NFκBIA may be important regulatory genes in calves' immune response to respiratory disease. Our findings provide a stable platform for a deeper

10
International Journal of Genomics understanding of the molecular mechanism behind calf BRD resistance in Xinjiang brown cattle.

Data Availability
The RNAseq data are available at NCBI GenBank database under accession number PRJNA702464.